On variance stabilisation by double Rao-Blackwellisation* 
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Abstract 

Population Monte Carlo has been introduced as a sequential importance sampling technique to 
overcome poor fit of the importance function. In this paper, we compare the performances of the 
original Population Monte Carlo algorithm with a modified version that eliminates the influence of 
the transition particle via a double Rao-Blackwellisation. This modification is shown to improve 
the exploration of the modes through an large simulation experiment on posterior distributions of 
mean mixtures of distributions. 

Keywords: Importance Sampling mixture, adaptive Monte Carlo, Population Monte Carlo, mul- 
timodal! ty, mixture of distributions, random walk. 
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1 Introduction 



When Cappe et al. (2004) introduced Population Monte Carlo (or PMC) as a repeated Sampling 



Importance Resampling procedure, the purpose was to overcome the shortcomings of Importance 
Sampling (abbreviated to IS in the following) procedures, while preserving the advantages of IS 
compared to alternatives such as Markov Chain Monte Carlo methods, namely the possibility of 
developing parallel implementations, which becomes more and more important with the generalisation 
of multiple core machines and computer clusters, and of allowing for easy assessment of the Monte 
Carlo error and, correlatively, for the development of on-line calibration mechanisms. 

Indeed, adaptive Monte Carlo naturally answers the difficulties of finding a "good" importance 
function in IS by gradually improving this importance function against a given target density, based 
on some form of Monte Carlo approximation. Previous simulations are used (see, e.g., Douc et al. 



2007a|b ) to modify proposal distributions represented as mixtures of fixed proposals in order to 
increase the weights of the most appropriate components. When the proposal is inspired from random 
walk structures as in Metropolis-Hastings algorithms and when the update of those weights is too 
crude, as in Cappe et al. (2004), the improvement is so short-sighted that multiple iterations do 



not increase the efficiency of the method, unless a Rao-Blackwell step is added to replace the actual 



proposal with the mixture proposal in the importance weight (see Douc et al. 2007a). More recently, 
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through the 2006-2008 project Adap'MC. The third author is particularly grateful to Andrew Barron for his suggestion 
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1 



Cappe et al. (2007) have thus replaced random walk proposals based on earlier samples with a new 
PMC scheme with parameterised proposals whose parameters are estimated from earlier samples, 
leading to a monotonic decrease in the entropic distance to the target distribution. There is indeed a 
myopic feature in random walk proposals, namely that, once a (first) sample has exhibited some high 
density regions for the target distribution, the algorithm is reluctant to allow for a wider exploration 
of the support of the target and it may thus miss important density regions, missing the energy to 
reach forward to these other regions of interest. 

In this paper, we nonetheless focus on random walk proposals because those kernels are more 
open to complex settings than on independent proposals — the later indeed require some preliminary 
knowledge about the target or else an major upgrade in computing power to face a much larger 
number of components in the mixture. More precisely, we present an experimental assessment of the 
use of a so-called double Rao-Blackwellisation scheme towards a better exploration of the modes of the 
target distribution. The second Rao-Blackwellisation step used in this double Rao-Blackwellisation 
scheme is essentially an integration over the particles of the previous sample used in the random 
walk proposal. While this leads to an additional computing cost in the derivation of the importance 
weights, double Rao-Blackwellisation undoubtedly brings a significant improvement in the stability 
of those importance weights and therefore justifies (in principle) the use of this additional step. 
Nevertheless, since an analytic proof of this improvement brought by double Rao-Blackwellisation 
is too delicate to contemplate, we use an intensive Monte Carlo experiment to establish the clear 
improvement in the case of a posterior distribution associated with a Gaussian location mixture and 
a sample with outliers — in order to increase the number of modes. 

The paper is organised as follows: In Section [2] we recall the basics of our population Monte Carlo 
algorithm and we introduce the double Rao-Blackwellisation modification. In Section [3j we motivate 
the choice of Gaussian mean mixtures as a worthy Monte Carlo experiment to test the mode finding 
abilities of both algorithm. Section [4] describes the implementation of the Monte Carlo experiment 
and describes the results. Section [5] contains a short discussion. 

2 Population Monte Carlo 

Given a target distribution with density n that is known up to a normalising constant, the grand 
scheme of PMC is the same as with other Monte Carlo — including MCMC — methods, namely to 
produce a sample that is distributed from tt without resorting to direct simulation from tt. Let us recall 



that, once a sample (X\, . . . ,Xn) is produced by sampling importance resampling (see, e.g., Robert 



and Casella 2004 Marin and Robert 2007), i.e., by first simulating X{ ~ f(%), second producing 
importance weights uii oc ir(Xi)/ f(Xi), and third resampling the Xi's by multinomial/bootstrap 
sampling with weights Ui, this SIR sample provides an approximation to the target distribution tt 
and can be used as a stepping stone towards a better approximation to tt. 

2.1 PMC basics 



In the PMC algorithm of |Cappe et TaL] ( |2004[ ) , if 

(X%, . . .,X N ) 

is a sample approximately distributed from tt, it is modified stochastically using an arbitrary Markov 
transition kernel q(x, x') so as to produce a new sample 



(X lf ... ,X r 



N ) 
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as X[ ~ q(Xi,x). Conducting a resampling step based on the IS weights Ui = Tr(X' i )/q(Xi, X-), we 
then produce a new sample 

(Xi, . . .,Xn) 

that equally constitutes an approximation to the target distribution %. It is however necessary to 
stress that, as established in Douc et al. (2007b), repeating this scheme in an iterative manner is 
only of interest if samples that have been previously simulated are used to update (or adapt) the 
kernel q(x, x'): a kernel q that is fixed over iterations does not modify the statistical properties of the 
samples. 



The choice made in Douc et al. (2007b) of an adaptive proposal kernel q represented as a mixture 
of fixed transition kernels qd, 



D D 
q a (x,x') = ^2a d q d (x,x') , ^a d =l, 
d=l d=l 



(1) 



can improve the efficiency of the kernel in terms of deviance (or relative entropy) from the target 
density. To achieve such an improvement, the weights a%, . . . , ud must be tuned adaptively at each 
iteration of the PMC algorithm. 



2.2 Weight updating 

If a 1 = fa* , . . . , a^j denote the mixture weights at the i-th iteration of the algorithm (where 
t = 1, . . . , T), the update of the cr's of Douc et al. (2007b) takes advantage of the latent variable 
structure that underlines the mixture representation, resulting in an integrated EM (Expectation- 
Maximisation) scheme. In the current setting, the latent variable Z is the standard component 
indicator in the mixture (see, e.g., Marin et al. 2005), with values in {1, . . . , D} such that the joint 
density / of x' and z given x satisfies 

f(z) = a z and f(x'\z, x) = q z (x, x') . 

The updating mechanism for the ct^'s then corresponds to setting the new parameter a t+1 equal to 

argmax E^f [Ef t {log(a z q z (X, X'))\X, X'}] = argmax e££' [e£ {log(a z )\X,X'}] , 

where the inner expectation is computed under the conditional distribution of Z for the current value 
of the parameter, a*, i.e. 

/ d 



with solution 



f(z\x, x', a*) = a* q z (x, x') / ^ a d q d (x, x') 

' d=\ 

a* +1 =E^' [f{d\X,X',a?)] . 



(2) 



2.3 Single Rao— Blackwellisation updates 

Given that ^ is not available in closed form, the adaptivity of the proposed procedure is achieved by 
approximating this actualising expression based on the sample that has been previously simulated. 
The implementation of the standard PMC algorithm relies on an arbitrary initial proposal ij,q that 
produces a pre-initial sample 

(-Xi,-ij • • • , -Xjv,-i) , 
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with importance weights ir(Xi-i) / fj,o(Xi _i). This initial sample allows for the derivation of a sample 



(Xi-i, . ■ . , X 



approximately distributed from n, using importance sampling based on those weights. The algorithm 
then picks arbitrary starting weights a d in the iV-dimensional simplex to produce a (truly) initial 
sample 

(Xifi, • • • , Xn,o) 



from the mixture 



D 



In the detail of its implementation, the production of this initial sample is naturally associated with 
latent variables (.Zi,o)l<i<iV that indicate from which component d of the mixture the corresponding 
Xifi has been generated [i = l,...,ra). From this stage, Douc et al. (2007b I proceed recursively. 
Starting at time t from a sample 

(Xij, • • • , XN,t) , 

associated with (Zij)i<i<N and with the current value of the weights a t,N , the normalised importance 
weights of the sample point Xi t is provided by 



Yld=i a d ,N( id(Xi it -i,Xi t t) 



■ N 

E 



7r(X 



i,t) 



(3) 



To approximate the update step ^ in practice, an initial possibility is 



N 



a t+i, = y^u^i^Zjj = d} 



i=l 

where the sum does not need renormalisation because the u>i t sum up to 1 (over i). 

This choice is however likely to be highly variable when N is small and/or the number of 
components D becomes larger. To robustify this update, Cappe et al. (2007) introduce a Rao- 
Blackwellisation step (see, e.g., |Robert and Casella 2004, Section 4.2) that consists in replacing the 
Bernoulli variable t{Zi t t = d} with its conditional expectation given Xi^ and Xi t—i, that is, 

TV 



a 



t+l,N 



E w m/ ( d 



Xi f . X, 



t,N 



i=l 



This replacement does not involve a significant increase in the computational cost of the algorithm, 
while both approximations are converging to the solution of Q as N grows to infinity. 



2.4 Double Rao-Blackwellisation 



While Cappe et al. (2007) showed through two experimental examples that the above Rao-Blackwel- 
lisation step does provide a significant improvement in the stability of the PMC weights (and, thus, 
in fine, a reduction in the number of iterations), there remains an extra-source of variation in the 
importance weights Q, namely the dependence of u^t on the previous sample point (or particle) 
Xi t t-i- Even though this dependence illustrates the fact that is indeed generated from the 



mixture 

D 



E 



a d ' q d (Xij-i,x) 



=i 
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and is thus providing a correct IS weight, it also shows a conditioning on the result of a (random) 
multinomial sampling in the previous iteration that led to the selection of Xij-i with probability 
Qi t—i- In other words, by de-conditioning, it is also the case that the sample point Xn has been 
generated from the (integrated) distribution 

N D 

^2°j,t-i ^2 a d Qd(Xj,t-i,x) , (4) 

j=\ d=l 

when averaging over all multinomial outputs. This double averaging over both the components of 
the mixture and the initial sample points is the reason why we call this representation double Rao- 
Blackwellisation. 

The appeal of using Q is that not only does the averaging over all possible sample points provide a 
most likely stabilisation of the weights, but it also eliminates a strange feature of the original approach, 
namely that two identical values Xij and Xj t t = Xij could have different importance weights simply 
because their conditioning sample values Xn-i and Xjt—i were different. 

Naturally, the replacement of the importance sampling distribution in ^ by Q has a cost of 
0(-/V 3 ) compared with the original 0(iV 2 ), but this is often negligible when compared with the cost 
of computing 7r(Xj it ). (We also stress that some time-saving steps could be taken in order to avoid 
computing all the qd(Xj,t-i,Xi t t) by considering first the distance between X/,t-i and X% t t and keeping 
only close neighbours within the sum although the increase in computing time in our case did not 
justify the filtering.) In the following experiment, the additional cost resulting from the double Rao- 
Blackwellisation does not induce a considerable upsurge in the computing time, even though it is not 
negligible. We indeed found an increase in the order of three to five times the original computing time 
for the same number N of sampling values. This obviously fails to account for the faster stabilisation 
of the IS approximation resulting from using double Rao-Blackwellisation. Note also that double 
Rao-Blackwellisation does not remove the need to sample the Xj f—i's: indeed, when simulating each 
new Xi t t from Q, we need to first select which Xjj-i is going to be conditioned upon, then second 
determine which component is to be used. 



3 The Gaussian mean mixture benchmark 



As benchmark for our Monte Carlo experiment, we consider the case of a one-dimensional Gaussian 
mean mixture distribution, namely 



, x n ~ pj\f(/J,i,af) + (1 - p)Af(fJ,2, 0\) 



(5) 



with p, a\ and a\ known and 9 = (/^i, fj^) being the parameter of the model, as in Marin et al. (2005 1 . 
Using a flat prior on 8 within a square region, we are thus interested in simulating from the posterior 
distribution associated with a given sample {x\, . . . ,x n ). The appeal of this example is that it is 
sufficiently simple to allow for an explicit characterisation of the attractive points for the adaptive 
procedure, being of dimension two, while still illustrating the variety of situations found in more 



realistic applications. In particular as already explained in Marin et al. (2005), the model contains 



at least one attractive point that does not correspond to the global minimum of the entropy criterion 
as well as some regions of attraction that can eventually lead to a failure of the algorithm when the 
data (x\, . . . ,x n ) is not simulated from the above mixture but by clumps that induce more modes 
in the posterior. Figure [T] illustrates quite well the diversity of posterior features when changing the 
repartition of the sample (x\, . . . , x n ). We stress the fact that those samples, called artificial samples, 
are not resulting from simulations from ([5]) but from simulations from a normal mixture with five 
equal and well-separated components centred in \x\ = and ±H2, ±2^2) and with variances 0.1. This 
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Figure 1: Series of posterior distributions associated with (|5| represented as level sets for various 
artificial samples (cci, . . . ,x n ) and different values of n, m, fj>2, a\ and a\. 

choice was made in order to increase the potential number of modes on the posterior surface when 
modelling the data with §5§ . 

4 Monte Carlo experiment 

Given the target distribution defined by (jHl), we want to study the improvement brought by the 
double Rao-Blackwellisation Q in terms of mode degeneracy, as well as to ascertain the additional 
cost of using double Rao-Blackwellisation. We therefore study the performance of both single and 
double Rao-Blackwellisation PMC over a whole range of samples, with various values of n and of the 
parameters p, fj,2 and 02 (since we can always use fii = and a\ = 1 without loss of generality). 

4.1 Automatic mode finding and classification 

A first difficulty, when building this experiment, is to determine the number of modes on the posterior 
surface for the data at hand. We achieve this goal by first discretising the parameter space in (/ii, H2), 
and then recovering the modal points by a brute-force search for maxima and saddle-points over the 
grid and then identifying their basins of attraction. (The algorithm is based on multidirectional 
calls to the function turnpoints () of the pastecs package.) This is obviously prone to overlook 
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some small local modes but it also allows for a quick identification of the modal regions and of their 
exploration by the PMC algorithm. 

Table [T] illustrates the results of this procedure for a range of values of (n,p, ^2,^2) ■ When 
conditioning upon {^2^2) the number of modes is increasing in \ii (since the simulated samples 
include five clusters that are better separated) and slightly decreasing in 02 (for apparently the same 
reason). A similar table varying upon the pair (n,p) does not show much variation in the number of 
modes (which does make sense, since only the relative magnitudes of the modes are changing). 



(0-2,^2) 


1.0 


1.5 


2.0 


2.5 


3.0 


3.5 


4.0 


4.5 


5.0 


1.0 


2.053 


2.184 


2.450 


2.863 


3.265 


3.524 


3.641 


3.707 


3.742 




(0477) 


(0.486) 


(1.273) 


(1.085) 


(0.895) 


(0.754) 


(0.706) 


(0.668) 


(0.627) 


1.5 


2.037 


2.124 


2.135 


2.165 


2.294 


2.494 


2.731 


2.958 


3.182 




(0.707) 


(0.362) 


(0.404) 


(0.565) 


(0.580) 


(0.668) 


(0.696) 


(0.729) 


(0.706) 


2.0 


2.026 


2.067 


2.121 


2.059 


2.083 


2.132 


2.259 


2.546 


2.807 




(0.862) 


(0.337) 


(0.359) 


(0.342) 


(0.350) 


(0.369) 


(0.505) 


(0.648) 


(0.636) 


2.5 


1.976 


2.024 


2.107 


2.067 


2.073 


2.184 


2.472 


2.757 


2.893 




(1.027) 


(0.354) 


(0.344) 


(0.302) 


(0.313) 


(O.434) 


(0.570) 


(0.533) 


(0.500) 


3.0 


1.993 


1.982 


2.073 


2.203 


2.137 


2.361 


2.743 


2.933 


3.040 




(1.1 47) 


(0.359) 


(0.329) 


(0.599) 


(0.469) 


(0.580) 


(0.572) 


(0.518) 


(0.486) 


3.5 


1.979 


1.964 


2.076 


2.186 


2.367 


2.662 


2.973 


3.181 


3.354 




(1.315) 


(0.432) 


(0.415) 


(0.519) 


(0.938) 


(0.844) 


(0.712) 


(0.678) 


(0.719) 


4.0 


1.961 


1.939 


2.078 


2.209 


2.412 


2.871 


3.214 


3.466 


3.752 




(1.469) 


(0.507) 


(0.490) 


(0.500) 


(0.851) 


(0.922) 


(0.807) 


(0.806) 


(0.860) 


4.5 


1.858 


1.893 


2.046 


2.220 


2.524 


3.101 


3.418 


3.777 


4.073 




(1-471) 


(0.504) 


(0-445) 


( 0.540) 


(0.832) 


(1.034) 


(0.898) 


(0.910) 


(0.901) 


5.0 


1.698 


1.886 


2.063 


2.273 


2.678 


3.187 


3.673 


4.046 


4.295 




( 1.247) 


(0.605) 


(0.566) 


(0.612) 


(1.128) 


(1.013) 


(1.033) 


(0.970) 


(0.867) 



Table 1: Average number (and standard deviation) of identified modes on a collection of 1470 
samples. 



The code for both implementing both versions of PMC and for evaluating their mode-finding 



abilities was written in R (R Development Core Team 2006) and is available as 



http:\\www.ceremade.dauphine.fr\~xian\2RB.R. Further documentation is available as 
http : \\www . ceremade . dauphine . f r\~xian\2RB . R . pdf . 



4.2 Output 

The experiment was run on 4860 samples on seven different machines for a total of 238140 simulations, 
using the same machine for both single and double Rao-Blackwellised versions at a given value of 
the parameters in order to keep the CPU comparison sensible. As seen from Table |2j the CPU time 
required to implement the double Rao-Blackwellised scheme is five to three time higher than the 
original PMC scheme, the additional time decreasing with n. (This CPU time does not include the 
mode finding and storing steps that are required for the comparison. It only corresponds to the 
execution of a regular PMC run with 10 iterations.) Note also that the increase in computing time 
is certainly less than linear. Both single and double Rao-Blackwellisation PMC samplers were used 
on the same samples, with different values of p, 02 and ji2- 

Turning now to the central tables, Tables [3}j6j we can see that the influence of 02 on the de- 
tection of the modes is relatively limited, in opposition to the influence of fi2- As p,2 increases 
(recall that fi\ = 0), a larger fraction of modes gets undetected. Both after 5 and 10 iterations of 
PMC, the performances of the double Rao-Blackwellised scheme are superior to those of the single 



7 



n 




u Zri,D 


20 


3.60 


15.59 




( U. OZJ 


(l.lzj 


30 


3.65 


15.51 




( U. OZJ 


{l.ldj 


40 


3.68 


15.59 




( U. OZJ 


(1.1JJ 


50 


3 72 


1 5 61 




( U. OZJ 


(I. Id) 


1 00 


3 8Q 


15 74 




( (/. Oo j 


(l.lbj 


500 


5 1 3 


1 7 1 Q 

J. l • J. iJ 




(0.62) 


(1.26) 


1000 


6.85 


18.92 






(1.42) 



Table 2: Overall mean CPU times and their standard error (in seconds) vs. sample size. Each 
value is obtained by averaging over 4860 generated samples, for a range of parameters (p = .1, . . . , .8, 
[12 = 1, • • • ,5, cr 2 = 1, . . . ,5). 



Rao-Blackwellisation in terms of mode detection, by a factor of about 5%. Running the PMC al- 
gorithm longer clearly has a downward impact on the number of detected modes, even though this 
impact is quite limited. It however points out the major issue that importance sampling algorithms 
like PMC are mostly unable to recover lost modes. Another interesting feature is that single Rao- 
Blackwellisation suffers more from the loss of modes between 5 and 10 iterations, compared with 
double Rao-Blackwellisation. As expected, the later is more robust because of the average over all 
past values in the PMC sample. Figure [2] presents the evolution of the capture rate for the single and 
the double Rao-Blackwellisations as a function of [ii after 5 and 10 iterations, obtained by averaging 
Tables |3H6] over c"2 . 



If we consider instead the output of this simulation experiment decomposed by the values of n 
and p, in Tables [7f|10| a main feature is the quick deterioration in the exploration of the modes 
as n increases. This is obviously to be expected given that the more observations, the steeper the 
slopes of the posterior surfaces. Thus, for small values of n, both types of Rao-Blackwellised PMC 
algorithms achieve high coverage rates, but larger values of n lead to poorer achievements. Once 
again, the robustness against the number of PMC iterations is superior in the case of the double 
Rao-Blackwellisation. Figure [3] presents the evolution of the capture rate for the single and the 
double Rao-Blackwellisations as a function of the sample size n after 5 and 10 iterations, obtained 
by averaging Tables [7 - 10 over p. 



Figure [4] illustrates the superior performances of the double Rao-Blackwellisation strategy, com- 
pared to the single one in terms of modes detection. Those three graphs plot side by side for a given 
sample the repartition of the PMC samples in the modal domains. The colour codes are associated 
with the five possible variances used in the 5 kernel proposal. In those three experiments, the per- 
sistence of the samples in all modes, as well as the concentration in the regions of importance, is 
noticeable. Both rows of Figure [4] interestingly contain a ridge structure in one of the modal basins. 
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2.0 


2.5 


3.0 


3.5 


4.0 


4.5 


5.0 


1.0 


0.743 


0.714 


0.631 


0.536 


0.452 


0.409 


0.399 


0.386 


0.387 




(0.260) 


(0.256) 


(0.254) 


(0.259) 


(0.232) 


(0.212) 


(0.216) 


(0.210) 


(0.209) 


1.5 


0.697 


0.697 


0.612 


0.573 


0.515 


0.472 


0.425 


0.386 


0.350 




(0.272) 


(0.255) 


(0.226) 


(0.208) 


(0.179) 


(0.175) 


(0.157) 


(0.145) 


(0.121) 


2.0 


0.653 


0.714 


0.614 


0.564 


0.552 


0.521 


0.495 


0.445 


0.397 




(0.273) 


(0.257) 


(0.226) 


(0.182) 


(0.173) 


(0.151) 


(0.147) 


(0.147) 


(0.131) 


2.5 


0.666 


0.718 


0.602 


0.563 


0.540 


0.509 


0.452 


0.394 


0.378 




(0.288) 


(0.256) 


(0.220) 


(0.181) 


(0.155) 


(0.147) 


(0.142) 


(0.114) 


(0.117) 


3.0 


0.684 


0.715 


0.605 


0.545 


0.526 


0.477 


0.402 


0.371 


0.355 




(0.302) 


(0.255) 


(0.219) 


(0.193) 


(0.156) 


(0.151) 


(0.119) 


(0.108) 


(0.094) 


3.5 


0.706 


0.718 


0.611 


0.538 


0.500 


0.434 


0.373 


0.345 


0.322 




(0.311) 


(0.257) 


(0.225) 


(0.179) 


(0.170) 


(0.152) 


(0.119) 


(0.103) 


(0.083) 


4.0 


0.725 


0.719 


0.610 


0.535 


0.484 


0.409 


0.347 


0.318 


0.291 




(0.314) 


(0.258) 


(0.225) 


(0.187) 


(0.155) 


(0.151) 


(0.108) 


(0.094) 


(0.084) 


4.5 


0.756 


0.736 


0.630 


0.536 


0.466 


0.379 


0.330 


0.293 


0.273 




(0.307) 


(0.259) 


(0.237) 


(0.192) 


(0.159) 


(0.137) 


(0.111) 


(0.086) 


(0.091) 


5.0 


0.781 


0.731 


0.627 


0.528 


0.462 


0.370 


0.312 


0.280 


0.253 




(0.293) 


(0.261) 


(0.238) 


(0.190) 


(0.179) 


(0.136) 


(0.111) 


(0.095) 


(0.079) 



Table 3: Average detection rate after 5 iterations of the single Rao-Blackwellised PMC algorithm 
a function of ai and /U2. The number of samples per entry is 1470, obtained for 7 values of n and 
values of p. 



(02,^2) 


1.0 


1.5 


2.0 


2.5 


3.0 


3.5 


4.0 


4.5 


5.0 


1.0 


0.794 


0.801 


0.725 


0.618 


0.528 


0.474 


0.463 


0.438 


0.440 




(0.258) 


(0.247) 


(0.266) 


(0.285) 


(0.265) 


(0.241) 


(0.238) 


(0.232) 


(0.229) 


1.5 


0.747 


0.811 


0.753 


0.700 


0.627 


0.558 


0.487 


0.439 


0.386 




(0.275) 


(0.240) 


(0.252) 


(0.257) 


(0.256) 


(0.248) 


(0.223) 


(0.216) 


(0.174) 


2.0 


0.679 


0.832 


0.763 


0.708 


0.669 


0.613 


0.566 


0.500 


0.443 




(0.279) 


(0.241) 


(0.255) 


(0.253) 




(0.232) 


(0.212) 


(0.206) 


(0.191) 


2.5 


0.679 


0.832 


0.756 


0.697 


0.648 


0.586 


0.514 


0.447 


0.416 




(0.291) 


( 0.240) 


(0.257) 


(0.251) 




(0.220) 


(0.212) 


(0.185) 


(0.160) 


3.0 


0.697 


0.826 


0.758 


0.672 


0.625 


0.555 


0.451 


0.411 


0.386 




(0.304) 


(0.242) 


(0.256) 


(0.263) 


(0.237) 


(0.226) 


(0.179) 


(0.159) 


(0.137) 


3.5 


0.712 


0.807 


0.745 


0.657 


0.591 


0.488 


0.424 


0.376 


0.355 




(0.312) 


(0.251) 


(0.260) 


(0.257) 


(0.244) 


(0.212) 


(0.178) 


(0.143) 


(0.131) 


4.0 


0.728 


0.788 


0.740 


0.646 


0.571 


0.457 


0.386 


0.350 


0.317 




(0.312) 


(0.258) 


(0.260) 


(0.257) 


(0.238) 


(0.204) 


(0.153) 


(0.138) 


(0.123) 


4.5 


0.758 


0.795 


0.738 


0.627 


0.540 


0.431 


0.363 


0.318 


0.298 




(0.307) 


(0.255) 


(0.262) 


(0.249) 


(0.229) 


(0.199) 


(0.146) 


(0.125) 


(0.121) 


5.0 


0.782 


0.783 


0.709 


0.617 


0.532 


0.413 


0.341 


0.305 


0.273 




(0.292) 


(0.259) 


(0.264) 


(0.253) 


(0.241) 


(0.182) 


(0.141) 


(0.124) 


(0.105) 



Table 4: Same legend as Table [3] for the double Rao-Blackwellised PMC algorithm. 
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(^2,^2) 


l.U 


1 K 

±.o 




z.o 


o.U 


O.O 


4.U 


4.0 


O.U 


1.0 


0.707 


0.641 


0.570 


0.482 


0.413 


0.376 


0.365 


0.356 


0.355 




( u.zoyj 


(0.244) 


( U.zzoJ 


( U.zzUJ 


( u.iyyj 


(U.lo 1 j 


/n 1 wy ) 
(U.lo 1 j 


/n 101) 
( U.lol J 


/n 1 q k ) 
(U.looJ 


1.5 


0.656 


0.599 


0.537 


0.516 


0.478 


0.443 


0.403 


0.371 


0.340 




( U.zozJ 


/n OOO) 

( U.zzzJ 


(U.l /UJ 


(0.14b) 


/n 1 fin) 

( u.izyj 


/ n 1 on) 

( u.izyj 


/ n 1 ok) 
(U.lzoJ 


/ n 1 1 n) 

(u.nyj 


/ n 1 no) 
( U.lUzJ 


2.0 


0.638 


0.622 


0.533 


0.518 


0.512 


0.491 


0.474 


0.427 


0.384 




( u.zoyj 


( U. zol J 


(U.lo 1) 


//i 1 1 n) 

(u.iiy) 


(U.lloJ 


/n nn i ) 

( u. uy4 ) 


fn 1 nQ) 
(U.lUoJ 


/n 1 1 k ) 
(U.lloJ 


/ n 1 n'y) 
( U.1U 1 J 


2.5 


0.656 


0.641 


0.531 


0.510 


0.507 


0.483 


0.436 


0.386 


0.366 




( U.zooJ 


( U.zo I J 


(U.l 00) 


(0.111) 


/ n inn) 

( U.1UUJ 


/ n 1 no) 
( U.lUzJ 


(0.111) 


/n nn k ) 

( u.uyoj 


/n nno) 

( u. uyzj 


3.0 


0.673 


0.660 


0.542 


0.497 


0.499 


0.456 


0.388 


0.359 


0.346 




/n ono) 
( U.oUzJ 


( U.z4oJ 


( U.looJ 


( U.IjUJ 


(U.llzJ 


/n 11/) 

(O.114) 


/n nnQ) 

( u. uyoj 


/n nQ i ) 
( U.U04) 


/n n Qn ) 
( U. UoUJ 


3.5 


0.700 


0.660 


0.554 


0.499 


0.472 


0.416 


0.362 


0.335 


0.315 




//I Q1 1 ) 
(U.ollJ 


/ 'n oik) 

(U.Z4OJ 


( U. 1 ozj 


(U. Iz 1 J 


//I 1 <D<D) 
( U. 1 zzj 


//) 100) 
( U. 1 ZZJ 


( u. uyo j 


(u.uooj 


(U.U 1 1 J 


4.0 


0.723 


0.670 


0.560 


0.498 


0.462 


0.390 


0.337 


0.308 


0.284 




(0.312) 


(0.248) 


(0.186) 


(0.139) 


(0.126) 


(0.121) 


(0.095) 


(0.075) 


(0.073) 


4.5 


0.755 


0.700 


0.572 


0.496 


0.445 


0.365 


0.319 


0.285 


0.267 




(0.308) 


(0.255) 


(0.201) 


(0.142) 


(0.129) 


(0.126) 


(0.090) 


(0.076) 


(0.079) 


5.0 


0.780 


0.700 


0.578 


0.501 


0.438 


0.357 


0.303 


0.270 


0.249 




(0.293) 


(0.258) 


(0.208) 


(0.159) 


(0.149) 


(0.120) 


(0.101) 


(0.078) 


(0.068) 



Table 5: Average detection rate after 10 iterations of the single Rao-Blackwellised PMC algorithm 
as a function of 02 and [ii- 



(0-2,^2) 


1.0 


1.5 


2.0 


2.5 


3.0 


3.5 


4.0 


4.5 


5.0 


1.0 


0.786 


0.785 


0.712 


0.610 


0.521 


0.467 


0.453 


0.431 


0.432 




(0.259) 


(0.251) 


(0.268) 


(0.284) 


(0.262) 


(0.240) 


(0.236) 


(0.228) 


(0.227) 


1.5 


0.745 


0.795 


0.733 


0.679 


0.612 


0.543 


0.479 


0.431 


0.379 




(0.275) 


(0.246) 


(0.253) 


(0.255) 


(0.249) 


(0.240) 


(0.216) 


(0.208) 


(0.167) 


2.0 


0.677 


0.813 


0.736 


0.684 


0.647 


0.600 


0.555 


0.494 


0.437 




(0.280) 


(0.246) 


(0.256) 


(0.247) 


(0.242) 


(0.223) 


(0.207) 


(0.200) 


(0.188) 


2.5 


0.675 


0.809 


0.726 


0.665 


0.623 


0.567 


0.500 


0.435 


0.409 




(0.291) 


(0.247) 


(0.256) 


(0.244) 


(0.229) 


(0.208) 


(0.199) 


(0.169) 


(0.152) 


3.0 


0.694 


0.804 


0.723 


0.643 


0.606 


0.537 


0.443 


0.404 


0.380 




(0.304) 


(0.251) 


(0.257) 


(0.257) 


(0.227) 


(0.211) 


(0.167) 


(0.154) 


(0.134) 


3.5 


0.711 


0.783 


0.706 


0.624 


0.568 


0.477 


0.413 


0.372 


0.350 




(0.311) 


(0.257) 


(0.256) 


(0.244) 


(0.233) 


(0.202) 


(0.169) 


(0.135) 


(0.124) 


4.0 


Q.727 


0.764 


0.704 


0.609 


0.542 


0.446 


0.377 


0.346 


0.312 




(0.313) 


(0.259) 


(0.255) 


(0.241) 


(0.218) 


(0.193) 


(0.145) 


(0.131) 


(0.116) 


4.5 


0.757 


0.766 


0.693 


0.596 


0.518 


0.416 


0.355 


0.313 


0.293 




(0.308) 


(0.259) 


(0.258) 


(0.237) 


(0.215) 


(0.185) 


(0.139) 


(0.117) 


(0.116) 


5.0 


0.782 


0.764 


0.665 


0.589 


0.512 


0.404 


0.334 


0.300 


0.268 




(0.293) 


(0.260) 


(0.255) 


(0.238) 


(0.228) 


(0.174) 


(0.134) 


(0.120) 


(0.098) 



Table 6: Same legend as Table [5] for the double Rao-Blackwellised PMC algorithm. 
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Modes detection race \s. [I , 




12 3 4 5 



Figure 2: Evolution of the capture rate for the single and the double Rao-Blackwellised versions as 
a function of \i2 after 5 and 10 iterations, obtained by averaging Tables [3]-[6] over 02 • 



(P> n ) 


20 


30 


40 


50 


100 


500 


1000 


0.10 


0.602 


0.565 


0.542 


0.519 


0.478 


0.423 


0.411 




(0.288) 


(0.276) 


(0.274) 


(0.260) 


(0.243) 


(0.198) 


(0.184) 


0.20 


0.580 


0.551 


0.527 


0.519 


0.480 


0.431 


0.414 




(0.278) 


(0.263) 


(0.249) 


(0.246) 


(0.222) 


(0.193) 


(0.173) 


0.30 


0.583 


0.552 


0.527 


0.529 


0.491 


0.445 


0.424 




(0.276) 


(0.262) 


(0.243) 


(0.246) 


(0.219) 


(0.193) 


(0.167) 


0.40 


0.582 


0.558 


0.538 


0.530 


0.496 


0.455 


0.432 




(0.272) 


(0.259) 


(0.253) 


(0.241) 


(0.221) 


(0.200) 


(0.174) 


0.50 


0.602 


0.572 


0.564 


0.553 


0.533 


0.476 


0.455 




(0.283) 


(0.263) 


(0.256) 


(0.247) 


(0.230) 


(0.192) 


(0.172) 


0.60 


0.597 


0.572 


0.544 


0.533 


0.499 


0.442 


0.427 




(0.276) 


(0.266) 


(0.256) 


(0.242) 


(0.220) 


(0.174) 


(0.153) 



Table 7: Average detection rate after 5 PMC iterations of the single Rao-Blackwellised algorithm as 
a function of n and p. The number of samples is 2835, obtained for 7 values of fj,% and 5 values of 02 • 



(P,n) 


20 


30 


40 


50 


100 


500 


1000 


0.10 


0.652 


0.619 


0.601 


0.577 


0.542 


0.463 


0.445 




(0.291) 


(0.286) 


(0.292) 


(0.283) 


(0.283) 


(0.249) 


(0.232) 


0.20 


0.644 


0.629 


0.613 


0.604 


0.566 


0.482 


0.453 




(0.286) 


(0.284) 


(0.283) 


(0.283) 


(0.279) 


(0.250) 


(0.229) 


0.30 


0.653 


0.645 


0.630 


0.628 


0.588 


0.500 


0.467 




(0.288) 


(0.288) 


(0.285) 


(0.284) 


(0.280) 


(0.256) 


(0.228) 


0.40 


0.664 


0.662 


0.654 


0.643 


0.598 


0.513 


0.472 




(0.291) 


(0.289) 


(0.290) 


(0.281) 


(0.283) 


(0.260) 


(0.227) 


0.50 


0.668 


0.651 


0.648 


0.639 


0.614 


0.524 


0.482 




(0.298) 


(0.290) 


(0.290) 


(0.284) 


(0.275) 


(0.241) 


(0.211) 


0.60 


0.665 


0.651 


0.627 


0.623 


0.582 


0.482 


0.450 




(0.291) 


(0.288) 


(0.285) 


(0.276) 


(0.272) 


(0.222) 


(0.189) 



Table 8: Same legend as Table [7] in the double Rao-Blackwellised case. 
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(p> n ) 


20 


30 


40 


50 


100 


500 


1000 


0.10 


0.577 


0.541 


0.519 


0.495 


0.457 


0.408 


0.398 




(0.281) 


(0.266) 


(0.259) 


(0.243) 


(0.223) 


(0.175) 


(0.161) 


0.20 


0.545 


0.514 


0.492 


0.485 


0.447 


0.413 


0.402 




(0.265) 


(0.242) 


(0.222) 


(0.218) 


(0.183) 


(0.164) 


(0.152) 


0.30 


0.541 


0.510 


0.489 


0.480 


0.456 


0.424 


0.414 




(0.257) 


(0.233) 


(0.214) 


(0.201) 


(0.180) 


(0.159) 


(0.146) 


0.40 


0.536 


0.516 


0.496 


0.487 


0.458 


0.431 


0.418 




(0.252) 


(0.234) 


(0.219) 


(0.207) 


(0.178) 


(0.162) 


(0.143) 


0.50 


0.564 


0.542 


0.528 


0.517 


0.497 


0.459 


0.443 




(0.269) 


(0.246) 


(0.232) 


(0.219) 


(0.197) 


(0.166) 


(0.156) 


0.60 


0.544 


0.519 


0.499 


0.487 


0.463 


0.426 


0.420 




(0.258) 


(0.238) 


(0.224) 


(0.206) 


(0.180) 


(0.146) 


(0.134) 



Table 9: Average detection rate after 10 PMC iterations for the single Rao-Blackwellised algorithm 
as a function of n and p. 



(P; n ) 


20 


30 


40 


50 


100 


500 


1000 


0.10 


0.639 


0.608 


0.590 


0.564 


0.529 


0.458 


0.441 




(0.291) 


(0.285) 


(0.288) 


(0.277) 


(0.275) 


(0.243) 


(0.227) 


0.20 


0.622 


0.612 


0.594 


0.585 


0.549 


0.475 


0.448 




(0.284) 


(0.281) 


(0.278) 


(0.276) 


(0.267) 


(0.246) 


(0.220) 


0.30 


0.629 


0.623 


0.608 


0.604 


0.571 


0.493 


0.461 




(0.284) 


(0.283) 


(0.279) 


(0.277) 


(0.274) 


(0.250) 


(0.223) 


0.40 


0.641 


0.636 


0.628 


0.620 


0.582 


0.506 


0.468 




(0.287) 


(0.285) 


(0.286) 


(0.276) 


(0.274) 


(0.254) 


(0.224) 


0.50 


0.642 


0.630 


0.627 


0.619 


0.600 


0.518 


0.480 




(0.293) 


(0.284) 


(0.283) 


(0.277) 


(0.269) 


(0.237) 


(0.208) 


0.60 


0.641 


0.628 


0.608 


0.604 


0.569 


0.475 


0.447 




(0.287) 


(0.281) 


(0.278) 


(0.271) 


(0.264) 


(0.215) 


(0.184) 



Table 10: Same legend as Table [9] for the double Rao-Blackwellised PMC algorithm. 



VI .il'-s (IlilcIu'mi rill l" vs. n 
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Figure 3: Evolution of the capture rate for the single and the double Rao-Blackwellisations as a 
function of the sample size n after 5 and 10 iterations, obtained by averaging Tables [7- 10 over p. 
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Figure 4: Comparison between PMC samples using a single Rao-Blackwellisation (left) and a double 
Rao-Blackwellisation (right). 
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5 Conclusions 



The double Rao-Blackwellised algorithm provides a more robust framework for adapting general 
importance sampling densities represented as mixtures in the sense that the influence of the starting 
points diminishes when compared with the original algorithm. It is thus unnecessary to rely on a 
large value of the number T of iterations of the PMC algorithm: with large enough sample sizes TV" 
at each iteration — especially on the initial iteration that requires many points to counter-weight a 
potentially poor initial proposal — , it is quite uncommon to fail to spot a stabilisation of both the 
estimates and of the entropy criterion within a few iterations. 
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